% Fourier法解二次元扩散方程
% 参考：Biner 相场编程
% Gitee Repo

clc
clear


L=1;
dx=0.05;
dt=0.0001;

[x y] = meshgrid(-L:dx:L);
n=size(x,1);

D=1;
sigma = dt*D/(dx^2);

if mod(n,2) == 1
	k1 = 2*pi/(2*L)*((1:ceil(n/2))'-1);
	k2 = -flip(k1);
	k2(end,:)=[];
	kx = [k1;k2];
	clear k1
	clear k2
else
	error('i cannot handle it');
end
ky = kx;
k2 = kx.^2+ky'.^2;

u0 = zeros(n,n,1);
u1 = zeros(n,n,1);
fu0 = zeros(n,n,1);
fu1 = zeros(n,n,1);

u0 = exp(-20*((x-0.2*L).^2+y.^2));
fu0 = (fft2(u0));

for tick = 1:10000
    % 显式法
	fu1 = fu0 - dt*D*k2.*fu0;
    % 隐式法
    fu1 = fu0./(1+dt*D*k2);
    fu0=fu1;

    % 作为对比的有限差分法
    diff_i = u0(1:n-2,2:n-1)-2*u0(2:n-1,2:n-1)+u0(3:n,2:n-1);
    diff_j = u0(2:n-1,1:n-2)-2*u0(2:n-1,2:n-1)+u0(2:n-1,3:n);
    u1(2:n-1,2:n-1) = u0(2:n-1,2:n-1) + sigma*(diff_i+diff_j);
    u0 = u1;

	if mod(tick,100) == 0
		clf
		axis([-L  L -L L 0 1])
        axis equal
        caxis([0 1])
		hold on

		scatter3(x,y,real(ifft2(fu1)))
		mesh(x,y,u1)

        view([30 30])
		drawnow
		pause(0.5)
	end
end
